Open source FCI code for quantum dots and effective interactions 
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We describe OpenFCI, an open source implementation of the full configuration-interaction 
method (FCI) for two-dimensional quantum dots with optional use of effective renormalized in- 
teractions. The code is written in C++ and is available under the Gnu General Public License. The 
code and core libraries are well documented and structured in a way such that customizations and 
generalizations to other systems and numerical methods are easy tasks. As examples we provide a 
matrix element tabulation program and an implementation of a simple model from nuclear physics, 
in addition to the quantum dot application itself. 
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I. INTRODUCTION 

Quantum dots, nanometre-scale semiconductor devices 
confining a varying number of electrons, have been stud- 
ied intensely in the last two decades. Quantum dots 
are fabricated using essentially macroscopic tools, for 
example etching techniques, but the resulting confine- 
ment allows for quantum mechanical behaviour of the 
electrons. Many of the parameters are directly control- 
lable, thereby justifying the term "artificial atoms" or 
"designer atoms". These considerations explain the im- 
mense research activityon these systems. For a general 
introduction, see Ref. [I[ and references therein. 

A very common model is that of a parabolic quantum 
dot, in which N electrons are confined in an isotropic 
harmonic oscillator potential in d spatial dimensions, 
where d is determined by the semi-conductor environ- 
ment. Electronic structure calculations on the parabolic 
dot and similar systems are often carried out using the 
full configuration- interaction method (FCI), also called 
exact diagonalization [1]. The Hamiltonian is then 
projected onto a finite-dimensional subspace of the N- 
electron Hilbert space and diagonalized. Care is taken in 
order to exploit dynamical and discrete symmetries of the 
exact problem, such as conservation of angular momen- 
tum and total electron spin, in order to block-diagonalize 
the Hamiltonian matrix and reduce the computational 
complexity. 

In this article, we describe OpenFCI, a recently de- 
veloped open source C++ code implementing the FCI 
method for quantum dots 0). The code has a generic 
framework in the shape of library functions, thereby al- 
lowing easy customization and extension to other sys- 
tems and methods, e.g., three-dimensional quantum dots 
or the nuclear no-core shell model. 

OpenFCI implements a rcnormalization of the two- 
body interactions, a technique widely used in nuclear no- 
core shell model calculations. This allows for accelerated 
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convergence with respect to Slater determinant basis size 
0,0]. To the author's knowledge, no other available code 
provides such effective interactions for quantum dot sys- 
tems. The code can be easily modified to create effective 
interactions for almost many-body problem using a har- 
monic oscillator basis. 



The code is developed in a Linux environment using 
Gnu C++, and is readily portable to other environments 
and compilers. The Fortran 77 libraries Lapack and 
Arpack are required, but as these are available on a wide 
range of platforms, portability should not be affected. 
OpenFCI is released under the Gnu General Public Li- 
cense (Gnu GPL) H and is documented using Doxygen 
|6|. As an open source project, the code can freely be 
used and modified. 



The article is organized as follows: In Section HU the 
FCI method is introduced in the context of the parabolic 
quantum dot, where we also discuss the reduction of the 
Hamiltonian matrix by means of commuting operators 
and configurational state functions. In Scction HTTl we dis- 
cuss the effective two-body interaction. As the technique 
is likely to be unfamiliar to most readers outside the nu- 
clear physics community, this is done in some detail. In 
Section ITVl we discuss the organization and use of Open- 
FCI. We also give some results from example runs, and 
in particular an analytically solvable non-trivial model 
due to Johnson and Payne is considered Q, where the 
only modification of the parabolic quantum dot is the 
interaction. Finally, we conclude our article in Section 



Two appendices have been provided, Appendix [A] de- 
tailing the heavily-used centre-of-mass transformation 
and Appendix [Bl discussing the exact numerical solution 
of the two-electron quantum dot needed for the effective 
interaction scheme. 
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II. FCI METHOD 
A. Hamiltonian in occupation number formalism 

We consider N electrons trapped in an isotropic har- 
monic oscillator potential in d spatial dimensions. The 
electrons interact via the Coulomb potential given by 
U(rij) = X/rij, where = — rj\\ is the inter-particle 
distance and A is a constant. The quantum dot Hamil- 
tonian then reads 

N N 
i— 1 i<j 

where the second sum runs over all pairs 1 < i < j < N, 
and where Hq(i) is the one-body Hamiltonian defined by 

Ho(i) :=-\^ + \\\n\\ 2 . 
The interaction strength A is given by 

/ m* 1 e 2 

A = \/-fc3-^ ' 2 

V wrr e 47reo 

where e is the dielectric constant of the semiconductor 
bulk, e 2 /47reo ~ 1.440 eV-nm, and ui = h/m*a 2 , a being 
the trap size and length unit, and m* being the effective 
electron mass. Typical values for GaAs quantum dots 
are e — 12.3, m* = 0.067 electron masses, and a = 20 
nm, yielding A = 2.059. The energy unit is huj, in this 
case hui = 2.84 meV. 

Choosing a complete set {4>a(x)} a eA of single-particle 
orbitals (where x — (r, s) denotes both spatial and spin 
degrees of freedom, and a — (a, a) denotes both generic 
spatial quantum numbers a and spin projection quan- 
tum numbers a = ±1), H can be written in occupation 
number form as 

H = h b a i*a b , a + i U c b d a la a lr a d,ra c ^, 

a,b o abed err 

(3) 

where (a a ) creates (destroys) a particle in the orbital 
4> a (x). These operators obey the usual anti-commutation 
relations 

{a a ,a}p\ = 5a tj 3, {a ai af}} = 0. (4) 

For a review of second quantization and occupation num- 
ber formalism, see for example Ref. 8]. The single- 
particle orbitals are chosen on the form 

where {<p a ('r)} are spinless orbitals and x<t( s ) = &<j,s are 
spinor basis functions corresponding to the eigenstates of 
the spin-projection operator S z with eigenvalues cr/2. 

It is important, that since the single-particle orbitals 
{ip a (r)} ae A are denumerable, we may choose an ordering 



on the set A, such that A can in fact be identified with 
a range of integers, A ~ {0, 1, 2, • • • , L/2}. In most ab 
initio systems L is infinite, since the Hilbert space is 
infinite-dimensional. Similarly, a = (a, +1) is identified 
with even integers, and a = (a, — 1) with odd integers, 
creating an ordering of the single-particle orbitals <p a (x), 
and a is identified with an integer < 1(a) < L. 

The single-particle matrix elements and the two- 
particle elements it"J are defined by 

h a b := fo> a |J?ob&) = J ~^Jf)H ip b (r) d d r, 

and 

Kd ■= ('Pa<Pb\U(ri2)\¥c<fd) 

= A / (p a {ri)<Pb(r 2 )— <Pc(ri)<fd(r2) d d rid d r 2 (5) 
J ri2 

respectively. 

The spatial orbitals (p a (r) are usually chosen as eigen- 
functions of Hq, so that = 5 a .be a . 

The basis functions for iV-particle Hilbert space are 
Slater determinants \& ai ,a 2 ,--- ,«») defined by 

\$ ai ,... aN ) ==0^0^ ■■■al lN \-}, 

where |— ) is the zero-particle vacuum. In terms of single- 
particle orbitals, the spatial representation is 

1 N 

' pSSjv t=l 

where Sn is the group of permutations of N symbols. 
The Slater determinants are anti-symmetric with respect 
to permutations of both xi and on, so that the orbital 
numbers on must all be distinct to give a nonzero func- 
tion. Each orbital is then occupied by at most one parti- 
cle. Moreover, for a given set {c^}^ of orbitals, one can 
create Nl distinct Slater determinants that are linearly 
dependent. In order to remove this ambiguity, we choose 
only orbital numbers such that /(ccj) < I{otj) whenever 
i < j. 

It follows, that there is a natural one-to-one corre- 
spondence between Slater determinants with N parti- 
cles and integers b whose binary representations have N 
bits set. (If \A\ = L < oo, the integers are limited to 
< b < 2 L .) Each bit position k corresponds to an 
orbital <f> a (x) through k = 1(a), and the bit is set if 
the orbital is occupied. Creating and destroying parti- 
cles in )QW ) simply amounts to setting or clearing 
bits (possibly obtaining the zero-vector if a particle is 
destroyed or created twice in the same orbital), keeping 
track of the possible sign change arising from bringing 
the set {ai} on ordered form using Eqn. (J3]). Note that 
the vacuum |— ) corresponds to b = 0, which is not the 
zero vector, but the single state with zero particles. 



3 



R 



2 ~ 



= 7i = l n = 2 n = 1 n = 



= 1 n = 1 re = 



= n = 1 



. = n = 



H 1 h 



H — m 



-4-3-2-1 1 2 3 4 

FIG. 1: Structure of single-particle orbitals of the two- 
dimensional harmonic oscillator. Angular momentum and 
shell number/energy on axes, and nodal quantum number n 
at each orbital. Orbital n — 0, m — —3 in shell R = 3 is 
occupied by two electrons for illustration. 



B. Model spaces 

The FCI calculations are done in a finite-dimensional 
subspace V of the iV-particle Hilbert space, called the 
model space. The model space has a basis B of Slater de- 
terminants, and V has the orthogonal projector P given 
by 



P := 



(0) 



I*b)e8 



The configuration-interaction method in general now 
amounts to diagonalizing (in the sense of finding a few 
of the lowest eigenvalues of) the, in general, large and 
sparse matrix PHP. The only approximation we have 
made is the truncation of the iV-particle Hilbert space. 

The model space V is seen to be a function of the single 
particle orbitals tp a (f), whom we choose to be the eigen- 
functions of Hq, i.e., harmonic oscillator eigcnfunctions. 
These may be given on several equivalent forms, but it is 
convenient to utilize rotational symmetry of Hq to create 
eigenf unctions of the projection of the angular momen- 
tum L z . In d = 2 dimensions we obtain the Fock-Darwin 
orbitals defined in polar coordinates by 



<y3n,m(7, 0) 



1 



:e im V m| LH(r 2 )e- r /2 . 



(J) 



Here L k n (x) = {-l) n [n\/{n+ |m|)!] 1 / 2 L*(a;) is the nor- 
malized generalized Laguerre polynomial. The factor 
(— 1)™ is for convenience, see Appendix I A II The har- 
monic oscillator energy is 2n+ \ m\ + 1 and the eigenvalue 
of L z = —id I do is m. All eigenfunctions with the same 
energy 2n+\m\ + l =: R +1 span a single-particle shell. 
The single-particle orbitals are illustrated in Fig. [TJ 
For a Slater determinant ... we have 



N 



2 fl o(i)l*ai,...,a, 



l$ai,-,a ? 



with 



JV 



where Ri = 2rii + \m,i\, and 



where M = YliLi mi - 

To complete our definition of P , we let 



B = Br = I |$ax,-,a, 



N \ 
i=l J 



(8) 



where R is called the energy cut, for obvious reasons. As 
R — ► cxd, the whole Hilbert space is spanned, and the 
eigenpairs of PHP converge to those of H. 



C. Configurational state functions and block 
diagonality 

In order to reduce the complexity of the computa- 
tions, we need to exploit symmetries of H. First of all, 
[H, L z ] — 0, and it is obvious that also [H, S z ] — 0, where 
the spin projection operator S z is given by 



S z ■= \ ^^a\,a a a. 



The Slater determinants are eigenvectors of both L z and 
S z with eigenvalues M and s z = Yli=i a i/^i respectively. 
We obtain a natural splitting of the model space V into 
subspaccs with constant angular momentum M and spin 
projection s z , viz, 



M,s z ■ 



M.s, 



M s z 



The diagonalization of H can thus be done within each 
space Pm,s z separately, amounting to diagonalizing indi- 
vidual blocks Pm.s z HPm,Sz- 

The Hamiltonian ([3]) also commutes with total electron 
spin S 2 , [S 2 , S z ] — 0, given by 



S 2 



1 



{S+S- + SS+), 



with 



s± 



so that a common basis for S z and S 2 would lead to even 
smaller matrix blocks. 

The eigenvalues of S 2 are on the form s(s + 1), where 
< 2s < N is an odd (even) integer for odd (even) 
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N. For a joint eigenfunction of S z and S 2 , called a con- 
figurational state function (CSF), \s z \ < s. The Slater 
determinants are, however, not eigenfunctions of S 2 , but 
such can be constructed by taking linear combinations 
of a small number Slater determinants. For details on 
this algorithm, see Ref. Q. Suffice it to say here, that 
S 2 only couples Slater determinants with identical sets of 
doubly occupied orbitals (meaning that 4>( a ,+) ancl 0(a,-) 
are both occupied, as in Figure [J) and singly occupied 



orbitals (meaning that only one of < 



0,+) 



and 



are 



occupied) . It is easy to see that S 2 does not couple Slater 
determinants in Vm,s z to another Vm>,s' ■ Thus, we ob- 
tain the splitting 



In Appendix [A] we provide the details of the centre- 
of-mass transformation. One then obtains the follow- 
ing prescription for the interaction matrix elements w"^: 
Let a = (/xi,^i), b = (^2,^2), c = (^3,^3), and d = 
(/i4,^4) be the circular quantum number equivalents of 
the usual polar coordinate quantum numbers Hi and m^. 
Due to conservation of angular momentum, we assume 
mi + m,2 = TO3 + 7714; otherwise, the matrix element 
uf d = 0. Define M = ^1+^2, M' = ^3 + ^4, N = vi + v 2 , 
and N' = v% + 1/4. Since m"J is linear in A, we set A = 1 
without loss of generality. Now, 



M 

«S = E 

P=P0 



N 

rp(M)rp(M') T (N) T (N') r \p 
P,P-2 p',M4 / j o'.va n, 



p—q\ 

n+s 5 



(9) 



9=90 



We stress that all the mentioned operators commute with 
each other, viz, 

[H,Qi\ = [0,-,%] =0, 

with Hi E {L Zl S z , S 2 }. If a modified problem breaks, 
say, rotational symmetry, such that [H, L z ] 7^ 0, we may 
still split the model space into to the eigenspaces of S z 
and S 2 . 



D. Matrix elements of Coulomb interaction 



The remaining ingredient in the FCI method is the 
Coulomb matrix elements defined in Eqn. ([5]). These 
can be calculated by first expanding L^(x) in powers of 
x using 



E 

rn— 



(-1)' 



(n + k)\ 



(n — m)!(fc + m)!m! 



and evaluating the resulting integral term-by-term by an- 
alytical methods (Iol | . The resulting expression is a seven- 
fold nested sum, which can be quite time-consuming, es- 
pecially if a large number of Fock-Darwin orbitals occurs 
in the basis B. Moreover, the terms are fractions of fac- 
torials with alternating signs, which is a potential source 
of loss of numerical precision. 

We therefore opt for a more indirect approach, giv- 
ing a procedure applicable to a wide range of potentials 
U(r 12) in addition to the Coulomb potential. Moreover, 
it can be generalized to arbitrary spatial dimensions d. 
The approach is based on directly transforming the prod- 
uct functions <fia(ri)<pb(f2) to the centre-of-mass system, 
where the interaction U(r 12) only acts on the relative 
coordinate, and then transforming back to the lab sys- 
tem. This reduces the computational cost to a doubly 
nested sum, as well as the pre-computation of the centre- 
of-mass transformation and the relative coordinate inter- 
action matrix. Both can be done exactly using Gaussian 
quadrature. The transformations to and from the cen- 
tre of mass frame are unitary transformations, which are 
stable and will not magnify round-off errors. 



where n = min(p, q), s = M' — M , p' = p + M 1 — M , 
q' = q + N' - N. Moreover, p = max(M' - M, 0) and 
q =max(JV'-iV,0). 

Here, are centre-of-mass transformation coeffi- 

cients defined in Appendix [A] while the relative coor- 
dinate interaction matrix elements CjTjj/, n,n' > 0, are 
defined by 

C S : = (Vn,m(r,0)\U(V2r)\<p n , im (r,6)) 

= 2 / r 2 ^Ll^(r 2 )L l ™ l (r 2 )U(V2r)e- r2 rQK)) 
Jo 

Depending on U(ri2), the integral is best computed us- 
ing generalized half-range Hermite quadrature (see Ap- 
pendix [B] and Ref. [Ill]) or Gauss-Hermite quadrature. 
Weights and abscissa for quadratures are conveniently 
computed using the Golub-Welsch algorithm [l2T |. which 
only depends on the ability to compute the coefficients 
of the three-term recursion relation for the polynomial 
class in question, as well as diagonalizing a symmetric 
tri-diagonal matrix. 

Let p(r) be a polynomial, and let a < 2 and /3 be 
non-negative constants. Then 



U(r 12 ) = rf 2 P(»la)e~ /3r ? a 



(11) 



admit exact evaluations using generalized half-range 
Gauss-Hcrmitc quadrature. The Coulomb potential, 
Gaussian potentials, and the parabolic interaction 
— Xr 2 2 /2 of the analytically solvable model treated in 
Sec. IIV CI belong to this class of potentials. 

In the case of a = 1 and p(r) — q(r ) (i.e., an even 
polynomial), the integral is more convenient to eval- 
uate using standard Gauss-Hermite quadrature. The 
Coulomb interaction falls into this class. 

Of course, one may let p(r) be a non-polynomial func- 
tion as well and still obtain very good results, as long 
as p(r) is well approximated with a polynomial, e.g., is 
smooth. 
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III. EFFECTIVE INTERACTIONS 
A. Motivation 

The FCI calculations converge relatively slowly as 
function of the model space parameter R [3j, as the er- 
ror AE in the eigenvalue behaves like o(R~ k ) in general, 
where k — 0(1). This behaviour comes from the singular 
nature of the Coulomb interaction. 

In Ref . [|| , numerical results using an effective interac- 
tion were presented. This method is widely used in no- 
core shell model calculations in nuclear physics, where 
the nucleon-nucleon interaction is basically unknown but 
highly singular This so-called sub-cluster effective 
interaction scheme replaces the Coulomb interaction (or 
another interaction) U(rij) = X/rij with a renormalized 
interaction U(i,j) obtained by a unitary transformation 
of the two-body Hamiltonian that decouples the model 
space V and its complement [Hj]. Therefore, the two- 
body problem becomes exact in a finite number of har- 
monic oscillator shells. Loosely speaking, the effective 
interaction incorporates information about the interac- 
tion's action outside the model space. In general, U(i,j) 
is non-linear in A and not a local potential. 

Using the renormalized U(i,j), the many-body system 
does not become exact, of course, but U(i,j) will per- 
form better than the bare interaction in this setting as 
well. To the author's knowledge, there exists no rigorous 
mathematical treatment with respect to this, but it has 
nevertheless enjoyed great success in the nuclear physics 
community pi. Il4, Ha ], and our numerical experiments 
unambiguously demonstrate that the convergence of the 
FCI method is indeed improved drastically [3[ , especially 
for N < 4 particles. We stress that the cost of producing 
U(i,j) is very small compared to the remaining calcula- 
tions. 



B. Unitary transformation of two-body 
Hamiltonian 

We now describe the unitary transformation of the 
two-body Hamiltonian (i.e., Eqn. (fT]) or ([3]) with N = 2) 
that de-couples V and its complement. This approach 
dates back as far as 1929, when Van Vleck introduced 
such a generic unitary transformation to de-couple the 
model space to first order in the interaction [l|| [TtJ • 

Let P be given by Eqn. and let D = dim(7 3 ). The 
idea is to find a unitary transformation H = HZ of H 
such that 

(1 - P)HP = 0, 

i.e., H is block diagonal. This implies that H c g defined 
by 

H eS := PHP 



has eigenvalues identical to D of those of the full operator 
H. Since D is finite, i7 c ff is called an effective Hamilto- 
nian. 

Selecting Z is equivalent to selecting a set of effective 
eigenpairs {{Ekt l 1 ^! ))}&Li) where Ef. is an eigenvalue of 
H and {|^| )}£=! C V are the effective eigenvectors; an 
orthonormal basis for V . It is clear that Z is not unique, 
since there are many ways to pick D eigenvalues of H, 
and for each such selection any unitary D x D matrix 
would yield an eigenvector set. 

However, some choices are more natural than others, 
since the eigenvectors and eigenvalues are usually con- 
tinuous functions of A. We then select the D eigenvalues 
Ek(X) that develop adiabatically from A = 0. For the 
corresponding effective eigenvectors ^^(A)), we choose 
the orthonormal set that minimizes the distance to the 
exact eigenvectors {|^/s )}£=!, i.e., 

D 

{|*f»E=i == axgmin £|||* fc ) - |^>|| 2 , (12) 

where the minimization is taken over orthonormal sets 
only. The effective eigenvectors also turn out to be con- 
tinuous functions of A, so i7 e ff will also be continuous. 

Let U is the D x D matrix whose columns contain 
Pl^Pfe) in the chosen basis, and let V be the correspond- 
ing matrix containing l^^). Clearly, V is unitary, while 
U only approximately so. Equation (|12|) can then be 
written 

V := argmin trace[([/ - U'){U - U')\ (13) 

w 

where the minimum is taken over all unitary matrices. If 
U has singular value decomposition given by 

U = XYY\ (14) 

the solution V is given by 

V:=XY\ (15) 

If E = diag(2?i, • • • , Ed) is the diagonal matrix whose 
elements are the chosen eigenvalues, we have 

H cfi = VEVK 

See Ref. [l3[ for a thorough discussion of the above pre- 
scription for -f/eff- 

Having computed the two-body -f/eff, we define the ef- 
fective interaction [7(1,2) by 

2 

17(1,2) :=H eS -Pj2 H o(i)P, 

i=l 

which gives meaning solely in the model space. In second 
quantization, 

abed err 
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and the TV-body H e g becomes (cf. Eqn. ([T])) 



iV 



N 

E 



U(i,j), 



with occupation number formalism form (cf. Eqn. (J3])) 



7.6cd err 



(16) 

Now, i/eff is well-defined in the space of iV-body Slater 
determinants where no pairs of occupied orbitals con- 
stitute a two-body state outside the two-particle model 
space, since then the matrix element 5"J would be un- 
defined. A little thought shows us that if [7(1,2) was 
computed in a two-body energy cut space with parame- 
ter R, H c s is well-defined on the many-body model space 
with the same cut R. 



C. A comment concerning the choice of model 
space 

The two-body problem is classically integrable, i.e., 
there exists 2d — 1 constants of motion f2j, such that 
their quantum mechanical observables commute with H 
and each other, viz, 

[H, Hi] = [Oi,%] = 0, for all i, j. 

Indeed, the centre-of-mass harmonic oscillator He de- 
fined in Eqn. (jT5J) below and the corresponding centre-of- 
mass angular momentum provides two constants, while 
total angular momentum L z provides a third. 

Using the model space V defined by an energy cut, we 
have 



case is problematic, since in the limit A — > 0, the ex- 
act eigenfunctions that develop adiabatically are not all 
either in the model space or in the complement. The adi- 
abatic continuation of the eigenpairs starting out in V is 
thus not well-defined. 

We comment, that the model space V 1 is often used 
in both no-core shell model calculations and quantum 
dot calculations, but the effective interaction becomes, 
in fact, ill-behaved in this case. 



D. Solution of the two-body problem 

What remains for the effective interaction, is the com- 
putation of the exact eigenpairs {{Ek, \^k) )}&Li- We 
must also solve the problem of following eigenpairs adia- 
batically from A = 0. 

For the two-body Coulomb problem, analytical solu- 
tions are available only for very special values for A [HI]. 
These are useless for our purpose, so we must use numer- 
ical methods. 

A direct application of the FCI method using Fock- 
Darwin orbitals with a large R' > R will converge slowly, 
and there is no device in the method for following eigen- 
values adiabatically. As the eigenvalues may cross, select- 
ing, e.g., the lowest eigenvalues will not work in general. 

For the two-body problem, the Pauli principle leads to 
a symmetric spatial wave function for the singlet s = 
spin state, and an anti-symmetric wave function for the 
triplet s = 1 spin states. For the spatial part, we exploit 
the integrability of the system as follows. Define centre- 
of-mass coordinates by 

1 



R := T2 (n+r 2 ) 



and 



[P., tk] = o 

as well, which is equivalent [l3| to 

[ffeff, flz] - 0, 



(17) 



so that H e g is integrable as well. In particular, U(l, 2) is 
block-diagonal with respect to fij . 

If we consider the commonly encountered model space 
V defined by the Slater determinant basis B 1 given by 



max(i?i) < R} 



instead of Eqn. 



we will have 
[P',#c]^0, 



as is easily verified. Indeed, V' is not an invariant sub- 
space of the centre-of-mass transformation T defined in 
Appendix [3] Thus, [H' eS , He] ^ 0, so that the centre- 
of-mass energy no longer is a constant of motion! The 
symmetry-breaking of the effective Hamiltonian in this 



r:= ^(n-r 2 ). 

Using these coordinates, the two-body Hamiltonian be- 
comes 



H = H {R) + 
='■ He + H rc i 



H (r) + U(V2r;X) 



(18) 



where ri2 = v2r := v^lHI- We have introduced the 
parameter A explicitly in the potential in this equation. 
H is clearly separable, and the centre-of-mass coordinate 
Hamiltonian Hq is a trivial harmonic oscillator, while the 
relative coordinate Hamiltonian can be written as 



1 2 

—r 
2 



U(V2r;X), 



where in polar coordinates r = (r cos 9, r sin 9) we have 



ld_ d_ 

r dr dr 



d9 2 ' 
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Applying separation of variables again, the eigcnfunc- 
tions of H le i can be written 



where ti is trie nodcil qusrrturn number. Un,m 

(r) satisfies 

where 



m\ 



1 d d m 2 1 9 . /- , . . . 
2r Or Or 2r z 2 



Equation (fT5|) is an eigenvalue problem in the Hilbert 
space L 2 ([0, oo), rdr), where the measure rdr is induced 
by the polar coordinate transformation. Although it is 
natural to try and solve the radial problem using Fock- 
Darwin orbitals, this will converge slowly. The solution 
to this problem is to use a radial basis of generalized half- 
range Hermite functions [ll[ . In Appendix [B] this is laid 
out in some detail. 

Equation (| 1 9j> is a one-dimensional equation, so there 
will be no degeneracy in the eigenvalues (i m .n for fixed 
to. In particular, the eigenvalues as function of the inter- 
action strength A will not cross, and will be continuous 
functions of A. We thus have n m ,n < Mm,n+i f° r au n > 
where n is the nodal quantum number. 

At A = we regain the harmonic oscillator eigen- 
values 2n + \m\ + 1. Correspondingly, the eigenfunc- 
tions ip m ,n(r,9) approaches the Fock-Darwin orbitals 
L Pm,n{i' 1 0) 1 i.e., the harmonic oscillator eigenfunctions. 
For the radial part, 

limn m ,„(r) = ff M(r) := V2r^ (r^ 2 / 2 . 
x — o 

Reintroducing spin, the full eigenfunctions = 
*ni,ttu,n2,m 3 are on t ne form 



2tt 



where s — for odd to-2, and s — 1 for even TO2, and 
\s z \ < s is an integer. 

Let Ri = 2rii + \mi\ be the shell numbers for the centre- 
of-mass coordinate and relative coordinate, respectively. 
The eigenvalue E = E numunaima is 

E = Rl + 1 + jU n2 ,m 2 ) 

with limit 



obeying R\ + < R. Turning on the interaction adi- 
abatically, the eigenpairs we must choose for the effec- 
tive Hamiltonian at a given A are exactly those with 
i?i + i?2 5: R- 

The model-space projection P^B needed in Eqns. (fT2|) 
and l|13p is now given by 



P^(x 1 ,x 2 ) 
where 

,|m| 



<Pni,mx{R)- 



n=0 



PR—Rl Un 2 ,"12 ( r ) 



7? 



X-s 



(21) 



where [^J is the integer part of a;. This operator thus 
projects onto the n + 1 first radial basis functions with 
given \m\. 

Due to Eqn. (fl7|) . the unitary operator Z can be de- 
composed into its action on blocks defined by tuples of 
n\,m\ and TO2 [13j]. The minimization (|12[) can then 
be applied on block-per-block basis as well. Each sub- 
problem is equivalent to the calculation of an effective 
Hamiltonian K c g of the radial problem for a given TO2 
and n. 

To this end, let n and m = to 2 be given. Let U be the 
(n + 1) x [n + 1) matrix whose elements are given by 

£4,fc = (slT 1 \u m , n ), 0<n,k <n, 

i.e., the model space projections of the exact eigenvec- 
tors with the lowest eigenvalues. Let U = XTY* be the 
singular value decomposition, and let V = XY^ . Then, 

K eS =Vdi a g(E ,--- ,E n )V^ 



and 

£in,\m\ 



*T e fr-diag(H+l,2- 



,2n- 



■1) 



is the (ni, toi, TO2)-block of the effective interaction. If 
we return to Eqn. ((3]), the effective interaction matrix 
elements u"^ are now given by replacing the matrix ele- 
ments C, 



n.n+s by the matrix elements C' n n , . 



Sn,\p-q\ 



where 



R- Ri - ItoI 



Ri =N + M -(p + q), 



where N, M, p and q are defined immediately after 
Eqn. ©. 



IV. CODE ORGANIZATION AND USE 



A. Overview 



E — ► iZj + R 2 + 2, 

which is the harmonic oscillator eigenvalue. 

As the centre-of-mass coordinate transformation con- 
serves harmonic oscillator energy, at A = 0, the eigen- 
functions that are in the model space are exactly those 



The main program is called QDOT, and processes a 
textual configuration file with problem parameters be- 
fore proceeding with the diagonalization of the Hamilto- 
nian. Eventually, it writes the resulting data to a Mat- 
lab/Gnu Octave compatible script for further process- 
ing. 
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As a C++ library as well as stand-alone application, 
OpenFCI is organized in several namespaces, which 
logically separate independent units. There are three 
main namespaces: manybody, gauss, and quantumdot. 
Put simply, manybody provides generic tools for many- 
body calculations, such as occupation number formal- 
ism, Slater determinants and CSFs, while gauss provides 
tools for orthogonal polynomials and Gaussian quadra- 
ture. These namespaces are independent of each other, 
and are in no way dependent on the particular quantum 
dot model. On the other hand, quantumdot synthesizes 
elements from the two former into a quantum dot FCI 
library. In QDOT, the main work is thus processing of 
the configuration file. 

Two other namespaces are also defined, being 
simple_sparse and simple_dense, which are, respec- 
tively, simple implementations of sparse and dense ma- 
trices suitable for our needs. We will not go into details 
in the present article. 

It should be clear that extending and customizing 
QDOT is a relatively easy task. The application QDOT 
is provided as a tool with a minimum of functionality, 
and the interested will almost certainly desire to further 
develop this small application. 

In order to help with getting started on such tasks, 
some stand-alone demonstration applications are pro- 
vided, all based on the core classes and functions. These 
include an interaction matrix element tabulator tabu- 
late, and a simple program PAIRING for studying the 
well-known pairing Hamiltonian fl9j |, which we will not 
discuss further here. Finally, there is a small interac- 
tive console-based Slater determinant demonstration pro- 
gram SLATER_demo as well. These applications will also 
serve as indicators of the flexibility of OpenFCI. 

OpenFCI does not yet support parallel computation 
on clusters of computers, using for example the Mes- 
sage Passing Interface [2CJ. Future versions will almost 
certainly be parallelized, but the present version in fact 
competes with parallel implementations of the standard 
FCI method with respect to convergence due to the ef- 
fective interaction implemented, see Sec. IIV CI The sim- 
ple structure of OpenFCI also allows users with less re- 
sources to compile and run the code. 



B. Core functionality 

The manybody namespace currently contains four main 
classes: Slater, CsfMachine, NChooseKBitset, and 
MatrixMachine. These will probably form the backbone 
of any manybody computation with OpenFCI. 

The class Slater provides Slater determinants, cre- 
ation and annihilation operators, and so on. It is based 
on the standard template library's (STL) bit set class, 
which provides generic bit set manipulations. The class 
NChooseKBitset provides means for generating sets of k 
objects out of n possible represented as bit-patterns, i.e., 
bit patterns corresponding to Slater determinants in the 



basis B or £>'. This results in a STL vector<Slater> ob- 
ject, which represent Slater determinant bases in Open- 
FCI. 

The class CsfMachine is a tool for converting a basis of 
Slater determinants into a basis of configurational state 
functions. These are represented as vector<csf _block> 
objects, where csf _block is a struct containing a few 
CSFs associated with the same set of Slater determinants 

A CSF basis is again input for the class 
MatrixMachine, which is a template class, and 
generates a sparse matrix PAP of an operator A, where 
P projects onto the basis. It also handles bases of pure 
Slater determinants as they are trivially dealt with 
in the CSF framework. The template parameter to 
MatrixMachine is a class that should provide the matrix 
elements hp, u"g, etc, of the generic operator given by 

a/3 a0"/5 
i 1 \ " a/37 t t t 

Notice, that the indices are generic orbitals, and not as- 
sumed to be on the form (a, a) as in Eqn. ([3]). 

Currently, only one-, two-, and three-body operators 
are implemented. The reason is, that the matrix ele- 
ments are not computed by directly applying the sum 
of creation- and annihilation operators to Slater deter- 
minants, since this approach, however natural, is very 
inefficient. Instead, we apply Wick's theorem directly Q 
on the matrix elements known to be not identically zero. 

In the gauss namespace, several functions are de- 
fined which computes sequences of orthogonal polyno- 
mials via recurrence relations and weights and abscissa 
for Gaussian quadratures based on these. The lat- 
ter is done using the Golub-Welsh algorithm, which 
only depends on being able to compute the coefficients 
of the recurrence relation [j^ . The most important 
functions are perhaps computeLaguerrePolysO and 
computeGenHalf GaussHermite () , which computes a se- 
quence of generalized Laguerre polynomials evaluated at 
a given set of points and quadrature rules for generalized 
half-range Hermite functions, respectively. 

Finally, the quantumdot namespace defines classes and 
functions that combined define the quantum dot prob- 
lem. The class RadialPotential encapsulates poten- 
tials on the form (jlip . It also computes effective inter- 
action blocks (7™'l m l. The class QdotHilbertSpace pro- 
vides means for generating the bases B and B' , utilizing 
conservation of angular momentum, using a fast, custom 
made algorithm independent of NChooseKBitset. The 
class QdotFci sews everything together and is basically 
a complete solver for the FCI method with effective in- 
teractions. 
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# ** Simple configuration file for qdot ** 

# — model space parameters — 

A = 4 # number of particles 

M = # angular momentum 

S = # total spin * 2 

R = 15 # cut of model space 



check of the validity of the calculations. 

By uncommenting the lines following the definition of 
lambda, we override the default Coulomb interaction, 
and produce a configuration file for the analytically solv- 
able model given by Johnson and Payne pj, where the 
Coulomb interaction is replaced by the parabolic inter- 
action 



# — interaction parameters — 
lambda =1.0 # interaction strength 



U(r 12 ) 



# pot_p = -0.5 

# pot_p_iseven = yes 

# pot_alpha = 

# pot_beta = 



# uncomment 

# these lines 

# for Johnson 

# Payne model 



# — computational parameters — 
use_energy_cut = yes 

nev =50 # no . eigenpairs 

use_veff = no 
matlab_output = results. m 

FIG. 2: Simple configuration file for QDOT 



TABLE I: Some ground state eigenvalues produced by QDOT 
for N = 3, 4 electrons with A = 2. Both the bare and the 
effective interaction are used 





N = 3, M 


= 0, a=\ 


N = 4, M = s = 


R 


E 




E 




6 


9.02370 


8.96523 


13.98824 


13.88832 


10 


8.97698 


8.95555 


13.86113 


13.83280 


14 


8.96800 


8.95465 


13.84491 


13.82848 


18 


8.96411 


8.95444 


13.83923 


13.82761 


22 


8.96191 


8.95435 


13.83626 


13.82730 


26 


8.96049 


8.95430 






30 


8.95950 


8.95428 







C. Sample runs 

A basic configuration file for QDOT is shown in Fig. [5] 
Varying the parameters lambda, R, S and the number of 
particles A |23[, and running QDOT each time, we pro- 
duce a table of ground state energies, shown in Table 
|TJ By changing the parameter use_vef f we turn on and 
off the effective interaction. The corresponding effective 
interaction ground states are also shown in the table. 
Notice, that with the effective interaction we obtain the 
same precision as the bare interaction, but with much 
smaller model spaces. This indicates that OpenFCI can 
produce results that in fact compete with parallel im- 
plementations of the standard FCI method, even in its 
present serial form. 

In Table [III we compare the ground state energies re- 
ported in Ref. @ with the alternative model space V' to 
the corresponding values produced by QDOT, also using 
T". This Table also appears in Ref. Q, and serves as a 



If A is sufficiently small, all the eigenvalues of this model 
are on the form 



E jjk = l+j + (k + N- l)Vl - NX, j,fc>0. (22) 

Since the potential is smooth, the eigenfunctions are all 
smooth, implying exponential convergence with respect 
to R Q. We therefore expect very accurate eigenvalues 
even with moderate R. In Table IIIII we show the first 
eigenvalues along with the error computed for 2V = 4 
electrons with A = 1/8. The computations are done in 
the M = s = s z = model space with R = 10 and 
R = 15. Some duplicates exist, and they are included 
for illustration purposes. It is evident, that the eigen- 
values become very accurate with increasing R; a clear 
indication of the correctness of the implementation. 



V. CONCLUSION AND OUTLOOK 

We have presented OpenFCI, an open source full con- 
figuration interaction implementation for quantum dots 
and similar systems. OpenFCI also implements a renor- 
malized effective interaction widely used in nuclear no- 
core shell model calculations, and we demonstrated that 
such interactions are indeed useful in the quantum dot 
calculations as well. 

OpenFCI is easy to extend and adapt. Possible appli- 
cations are computations on systems with more general 
symmetry-breaking geometries and in d = 3 spatial di- 
mensions. Also, a generalization of the CSF part of the 
code to handle isobaric spin would allow us to handle 
nuclear systems. 

There is one more symmetry of the Hamiltonian H that 
can be exploited, namely that of conservation of centre- 
of-mass motion, which would further reduce the block 
sizes of the matrices. We exploited this symmetry for the 
effective interaction, but it is a fact that it is a symmetry 
for the full Hamiltonian as well. Using the energy cut 
model space V we may take care of this symmetry in a 
way similar to the CSF treatment 21 1. 

As mentioned, we have not parallelized the code at 
the time of writing, but it is not difficult to do so. A 
future version will almost certainly provide parallelized 
executables, for example using the Message Passing In- 
terface 12011. 
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TABLE II: Comparison of current code and Ref. 0, taken from Ref. 0| 











R = 5 




R = 6 




R = 7 




N 


A 


M 


2s 


Current 


Ref. 9 


Current 


Ref. 9 


Current 


Ref. 9 


2 


1 





o 


tj . V7 -L *J V_» _■ 




S 01 1 020 




3 00Q236 






2 


o 


o 


3 733598 


3.7338 


3 731 057 

iJ. 1 Ul \J' J I 


3.7312 


3 729324 


3.7295 






1 


2 


4.143592 


4.1437 


4.142946 


4.1431 


4.142581 


4.1427 


3 


2 


1 


1 


8.175035 


8.1755 


8.169913 




8.166708 


8.1671 




4 


1 


1 


11.04480 


11.046 


11.04338 




11.04254 


11.043 









3 


11.05428 


11.055 


11.05325 




11.05262 


11.053 


4 


6 








23.68944 


23.691 


23.65559 




23.64832 


23.650 






2 


4 


23.86769 


23.870 


23.80796 




23.80373 


23.805 


5 


2 





5 


21.15093 


21.15 


21.13414 


21.13 


21.12992 


21.13 




4 





5 


29.43528 


29.44 


29.30898 


29.31 


29.30251 


29.30 



TABLE III: Results from diagonalizing the Johnson and 
Payne model. Many digits are included due to comparison 
with exact results and high precision 



R= 10 



R= 15 



E 


AE 




E 


AE 




4.535550207816 


1.63 


10" 


5 


4.535533958447 


5.25 


10" 




5.950417930316 


6.70 


10" 


-4 


5.949751427847 


3.96 


10" 


6 


5.950417930316 


6.70 


10" 


-4 


5.949751427847 


3.96 


10" 


6 


5.950417930316 


6.70 


10" 


-4 


5.949751427847 


3.96 


10" 


6 


5.951592166603 


1.84 


10" 


-3 


5.949760599290 


1.31 


10" 


5 


6.243059891817 


4.19 


10" 


-4 


6.242642740293 


2.05 


10" 


6 


6.243059891817 


4.19 


10" 


-4 


6.242642740293 


2.05 


10" 


6 


6.535776573577 


2.43 


10" 


-4 


6.535534873729 


9.68 


10" 


7 


6.535776573577 


2.43 


10" 


-4 


6.535534873729 


9.68 


10" 


7 


6.535776573577 


2.43 


10" 


-4 


6.535534873729 


9.68 


10" 


7 


7.375904323762 


1.19 


10" 


-2 


7.364103882564 


1.43 


10" 


4 


7.375904323762 


1.19 


10" 


-2 


7.364103882564 


1.43 


10" 


4 


7.375904323762 


1.19 


10" 


-2 


7.364103882564 


1.43 


10" 


4 


7.375904323762 


1.19 


10" 


-2 


7.364103882564 


1.43 


10" 


4 


7.375904323762 


1.19 


10" 


-2 


7.364103882564 


1.43 


10" 


4 


7.393706556283 


2.97 


10" 


-2 


7.364440927813 


4.80 


10" 


4 


7.393706556283 


2.97 


10" 


-2 


7.364440927813 


4.80 


10" 


4 


7.393706556283 


2.97 


10" 


-2 


7.364440927813 


4.80 


10" 


4 


7.410720999386 


4.68 


10" 


-2 


7.364876152101 


9.15 


10" 


4 


7.665921446569 


9.07 


10" 


-3 


7.656945606956 


9.14 


10" 


5 



APPENDIX A: CENTRE OF MASS 
TRANSFORMATION 



Cartesian coordinates 



In this appendix, we derive the centre-of-mass (COM) 
transformation utilized in Eqn. ([9]) for the interaction 



1 cd- 



The one-dimensional harmonic oscillator (HO) Hamil- 
tonian (p 2 x + x 2 )/2 is easily diagonalized to yield eigen- 
functions on the form 

= {n\)- l l 2 A n Mx), (Al) 
where A x := (x — ip x )/y/2 is the raising operator in the 



^-coordinate, and where (f>o(x) = n^ 1 / 4 cxp(— x 2 /2). The 
eigenvalues are n + 1/2. 

Using separation of variables, the two-dimensional HO 
Hq in x\ and X2 is found to have eigenfunctions on the 
form ® nun2 {xi,x 2 ) := ^1(^1)^2(^2) and eigenvalues 
n\ + ri2 + 1. Define the raising operators A Xi := (xi — 
ip x . ,)/V2, so that 

1V,n 2 (£) := (n.W.y^A^A^oix), (A2) 
where the non-degenerate ground state is given by 



$0,0(2:1, £2) = —=e 

\/7T 



-(xf+xD/2 



(A3) 



Note that this can equally well describe two (distinguish- 
able and spinless) particles in one dimension. 

To this end, we introduce normalized COM frame co- 
ordinates by 



1 



1 1 

1 -1 



F 



(A4) 



The matrix F is symmetric and orthogonal, i.e., F T F = 
F 2 = 1, transforming a set of Cartesian coordinates into 
another. The operator Hq is invariant under this trans- 
formation, so the eigenfunctions have the same form with 
respect to these coordinates, viz, 

^,„ 2 (6,6) := ^(£1)^(6) 

= (m\n 2 \)-^ 2 A^A^' 0t0 , (A5) 

where A^ = (£j — ip'A/^/2 are the raising operators with 
respect to the COM coordinates, and p\ are the corre- 
sponding momentum components. 
Define the operator T by 



Tip(x 1 ,X2) := ^(£1,6) 
so that 



J, 



xi + x 2 Xi — x 2 
V2 ' V2 



. (A6) 



(A7) 
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Since T maps eigenfunctions in the two frames onto each 
other, T must be a unitary operator, and the invariance 
of Ho under the coordinate transformation is the same as 
[Hq,T] — 0, i.e., that energy is conserved. This in turn 
means that T is block diagonal with respect to each shell 
R = ni + ri2, viz, 

R 

n=0 
R 



(A8) 



where is the (R + 1) X (R + 1) transformation ma- 
trix within shell R. It is real, symmetric, and orthogo- 
nal. Numerically, the matrix elements are conveniently 
computed using two-dimensional Gauss- Hermitc quadra- 
ture of sufficiently high order, producing exact matrix 
elements. 

In a two-dimensional setting, the two-particle har- 
monic oscillator becomes a 4-dimensional oscillator. Let 
fi = (xi,yi), i — 1,2, be the particles' coordinates, and 
let a = (7711,712) and b = (7772,712) to compress the nota- 
tion a little. An eigenfunction is now on the form 



*o,6(n,r 2 ) 



$a(ri)$ 6 (r 2 ) 



= CA™} A™* A™ 0j0 (fi , r 2 ) ( A9) 



Vi ' 



where C — (mi!ni!m2!772!) _1//2 . 

The COM coordinate transformation now acts in the 
x and y directions separately, viz, F acts on Xi and yi 
to yield the COM coordinates and 7^: [£i,£2] T = 
F[xi : X2\ T and [/71 ,772]^" = F[yi,y2] T - The induced oper- 
ator T again conserves energy. Let M = mi + m 2 and 
N = ni + 77 2 . It is readily verifiable that the COM frame 
transformation becomes 



T^S M-m 2 ,N-n 2 ,rn 2 ,n 2 '■— ^ M-m 2 ,N -n 2 ,m 2 ,n 2 



M N 
p=0 <;=0 



(A10) 



Note that the shell number is i? = N + M, which is 
conserved by T. 



2. Centre of mass transformation for Fock-Darwin 
orbitals 

Consider a Fock-Darwin orbital tp n>m {r) in shell R = 
2n + I m I with energy R + 1. It is straightforward but 
somewhat tedious to show that these can be written in 
terms of co-called circular raising operators B + and i?_ 
[H defined by 



Letting fi = n + max(0, m) and v = n + max(0, —777) 
(which gives /j, v > 0) one obtains 



(A12) 



which should be compared with Eqn. (|A2|) . Moreover, 
R = fi + v and m = ji — v, giving energy and angular 
momentum, respectively. We comment that this is the 
reason for the non-standard factor (—1)™ in the normal- 
ization of the Fock-Darwin orbitals in Eqn. 0. 
Let a two-particle HO state be given by 



Ml ifi <M2, v 2 ■ — l fin 1 ,m 1 {' r l) l Pn 2 ,m 2 
Ml T)V\ r>M2 r>v 2 



(r 2 ), 



CB$B?_B%_B?_ * ,o(ri , ?2 )(A13) 



where /i; = 77^ + max(0, mi) and i>i — + max(0, —mi), 
and where C = (/iil^il/j^!^!) We will now prove 

that, in fact, when applying the centre-of-mass transfor- 
mation to Eqn. (|A13[) . we obtain an expression on the 
same form as Eqn. (|A10|) viz, 



M N 



/ j M2.P / t v 2 ,q 
p=0 q=Q 



M — n 2 ,N — v 2 ,ii 2! v 2 
M-p,N-q,p,q, (A14) 



where M := \i\ + /^2 and N := v\ + t^. 

To this end, we return to the raising operators A^ and 
.A^ , and express them in terms of A Xi and A Vi . By using 
Eqn. I|A4|) . we obtain 



F 



(A15) 



and similarly for A 7]i in terms of A Vi . In terms of the 
raising operators, the COM transformation becomes 



^i.ni,m 2 ,n 2 — C (A Xl -\- A X2 ) (A yi + A y2 ) 



x (A Xl — A X2 ) 2 (A Wl 



where 



_ ( 9 nx+n 2 +mi+m 2 1 , 



C= (2 



77l!77 2 ! 777l!777 2 !) 1//2 , 



*^16) 



(A17) 



and we have used Eqn. (|A9j) . but in the analogous COM 
case. Expanding the powers using the binomial formula 
(and the fact that the raising operators commute), we 
obtain a linear combination of the individual eigenfunc- 
tions, which must be identical to Eqn. (jAlOj) . 

Let the COM circular ladder operators be defined by 



(A18) 





1 


'1 i ' 




"V 






1 -i 




A, h _ 



Using Eqn. (|A15[) , we obtain that the circular raising 
operators transform in the same way as the Cartesian 
operators when going to the COM frame, i.e., 





1 


'1 i ' 




A x 


B_ 




1 -i 




Ay. 



(All) 



B' 



= F 



^2H 



(A19) 
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and similarly for B[_ and terms of Bi^. Using 
Eqn. (|A13p in the COM case, we obtain 

*™ 2 ,„ 2 = C (B 1+ + B 2+ r (Bi- + B 2 .p 

x (B 1+ - B 2+ r [B x _ - B 2 _) U2 *6&?0) 

with 

C= {2^ +Ul+ ^ +v ^ l W.^ 2 W-y 1/2 ■ (A21) 

Eqn. (|A20[) is on the same form as Eqn. (|A16p . Again, 
by expanding the powers using the binomial formula (and 
that the raising operators commute), we obtain a linear 
combination with coefficients identical to those of the ex- 
pansion of Eqn. (|Al"6)) . It then follows that Eqn. (fAl~4|) 
holds. 



APPENDIX B: NUMERICAL TREATMENT OF 
RADIAL PROBLEM 

We now briefly discuss the numerical method used for 
solving the radial problem (|19|) . i.e., the eigenvalue prob- 
lem for the operator K\ m \ defined in Eqn. ([20)) . This is an 
eigenproblem in the Hilbert space L 2 ([0, oo), rdr), where 
the measure rdr is induced by the polar coordinate trans- 
formation. The inner product on this space is thus given 
by 

(f\g) = / f(r)g(r)rdr. (Bl) 
Jo 

Let the Fock-Darwin orbitals be given by 

p imQ 

y ntm {r,9) = -j=gW{r), (B2) 

with radial part 

fl H(r) := ^4"V> |m| exp(-r 2 /2) (B3) 

Thus, 

and these functions form an orthonormal sequence in 
i 2 ([0, oo), rdr) for fixed \m\. 

In the electronic case, U(V2r) = X/V2r has a singular- 
ity at r = which gives rise to a cusp in the eigenfunction 
u mtn (r) at r = 0, or in one of its derivatives. Away from 
r = 0, the eigenfunction is smooth. These considera- 
tions are also true for more general potentials smooth for 
r > 0. 

Diagonalizing the matrix of K with respect to the 

truncated basis {5n m '(r)}™ =0 will give eigenpairs converg- 
ing slowly with respect to increasing n due to the non- 
smoothness of tin m (r) at r = 0. This is easily seen for 
the m = case and A = 1, which has the exact ground 
state 



uo,o( r ) = 




a polynomial of odd degree multiplied by a Gaussian. 
The cusp at r = is evident. However, 

which are all even polynomials. It is clear, that n must 
be large to resolve the cusp of Mo,o(f)- 

The eigenproblem is best solved using a basis of gener- 
alized half-range Hermite functions fj (r) [ll| , which will 
resolve the cusp nicely. These functions are defined by 

fj(r) :=P,(r)exp(-r 2 /2), 

where Pj(r) are the orthonormal polynomials defined 
by Gram-Schmidt orthogonalization of the monomials r k 
with respect to the weight function rexp(— r 2 ). Thus, 

/>oo 

</il/i'>= / f j (r)fAr)rdr = 8 jd , 
Jo 

The fundamental difference between fj(r) and (7« m '(r) 
is that the latter contains only even (odd) powers of r 
for even (odd) |m|. Both sets constitute orthonormal 
bases, but fj(r) will in general have better approximation 
properties. 

Moreover, since deg(rl m 'xll rl ' (r 2 )) = 2n + \m\, 

2n+\m\ 

at\r)= (B4) 
j=o 

gives the Fock-Darwin orbitals as a finite linear combi- 
nation of the generalized half-range Hermite functions, 
while the converse is not possible. 

Computing the matrix of K with respect to {/j(r)}j_ 
and diagonalizing will give eigenpairs converging expo- 
nentially fast with respect to increasing j. The resulting 
eigenf unctions' expansion in gl^ are readily computed 
using Eqn. (|B4[) . whose coefficients (fj\gn) can be com- 
puted numerically exactly using Gaussian quadrature in- 
duced by -Pj(r), for J sufficiently large. 

The basis size j to use in the diagonalization depends 
on how many eigenfunctions n we desire. We adjust j 
semi-empirically, noting that 2n + \m\ is sufficient to re- 
solve f/li" 1 ', and assuming that the exact eigenfunctions 
are dominated by the latter. We then add a fixed num- 
ber jo to get j — 2ft+|m|+jo, and numerical experiments 
confirm that this produces eigenvalues that indeed have 
converged within desired precision. 
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